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Abstract. We use an iterative generalized least squares map-making algorithm, in conjunction with Monte Carlo 
techniques, to obtain estimates of the angular power spectrum from cosmic microwave background (CMB) maps. 
This is achieved by characterizing and removing the instrumental noise contribution in multipole space. This 
technique produces unbiased estimates and can be applied to an arbitrary experiment. In this paper, we use it 
on realistic simulations of Planck Low Frequency Instrument (LFI) observations, showing that it can lead to fast 
and reliable estimation of the CMB angular power spectrum from megapixel maps. 
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1. Introduction 

State of the art measurements of the cosmic microwave 
background (CMB) anisotropy are affected by non negligi- 
ble and (in the case of one horned experiments) correlated 
instrumental noise. In order to minimize this contaminant 
one resorts to non trivial statistical techniques which al- 
most invariably require knowledge of the correlation struc- 
ture of underlying detector noise, to be measured from the 
data themselves. As a consequence, methods to estimate 
the noise correlation properties out of time ordered data 
(TOD) have been proposed (Dore et al. 2001; Stompor et 
al. 2002; Natoli et al. 2002). 

Most cosmological information is encoded in the an- 
gular power spectrum of the CMB anisotropics. However, 
the size of modern CMB datasets and the presence of cor- 
related instrumental noise in the observations make it un- 
feasible to extract the power spectrum from CMB maps 
by performing standard and robust matrix manipulations 
commonly used to solve linear systems. Computationally, 
obtaining a brute force maximum likelihood estimation of 
the power spectrum requires at least TVp operations, where 
Afp is the number of pixels in the map. Current CMB maps 
from balloon-borne experiments have M v ~ 10 4 — 10 5 . 
Brute force power spectrum estimation from these maps is 
already prohibitive, and will become totally unattainable 
for upcoming space missions such as NASA's MAP satel- 
litef] or ESA's Planck Surveyor^, whose maps will have 

K > io 6 . 
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A number of strategies have been proposed in the 
past to address the problem of power spectrum estima- 
tion in a computationally feasible way. The MADCAP 
package (Borrill 1999) uses a parallel implementation of 
a quadratic estimator technique (Bond et al. 1998) to 
obtain maximum likelihood estimates of the power spec- 
trum. This method, however, will certainly be too time 
consuming for future satellite data sets. Dore et al. (2001) 
applied the same approach in a hierarchical fashion on 
subsets of large data sets, lowering somewhat the compu- 
tational time requirements. Szapudi et al. (2001) adopted 
an entirely different strategy, extracting the power spec- 
trum from the 2-point correlation function of the map. 
Although this approach is applicable to a generic data set 
and might in principle account for noise correlations, until 
now it has only been tested in the case of uniform white 
noise. Finally, other methods were based on simplifying 
assumptions tailored on specific experimental strategies 
(e.g., Oh et al. (1999) for the MAP satellite; Wandelt k 
Hansen (2001) for the Planck Surveyor). 

In this work, we focus on characterizing and remov- 
ing the instrumental noise contribution in multipole space, 
in order to obtain unbiased estimates of the CMB power 
spectrum. This approach is based on using map-making 
techniques to project estimates of the time stream noise 
on the sky, according to the experimental scanning strat- 
egy, and on Monte Carlo (MC) simulations to determine 
the statistical properties of the noise in multipole space. 
Rather than being based on a noise model or an approxi- 
mation (such as, e.g., uncorrelated noise) the time stream 
noise properties are estimated directly from the data. 
A similar strategy was developed independently in the 
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MASTER package by Hivon et al. (2001) and applied to 
the analysis of the BOOMERanG data (Netterfield et al. 
2002). In this case, fast projection of the noise on the sky 
was achieved through non optimal map-making, and filter- 
ing of the data in time domain was required to reduce the 
correlated noise contribution. As such, MASTER does not 
attempt to create a useful map as an end-product. Instead, 
we use the iterative generalized least squares (IGLS) map- 
making algorithm described in Natoli et al. (2001) which is 
fast enough to allow for the generation of an appropriate 
number of simulations in a reasonable time. This algo- 
rithm has optimal statistical properties (i.e. produces un- 
biased and minimum variance estimates of the map) and 
hence it does not require any filtering of the time stream. 
As an application, we use realistic simulations of Planck 
Low Frequency Instrument (LFI) observations to show 
that we can successfully recover the CMB power spec- 
trum even from megapixel maps. This is, to date, the first 
computationally feasible and unbiased pixel-based power 
spectrum estimator for Planck which uses information on 
the correct (i.e. estimated from the data themselves) in- 
strumental noise covariance. 

This paper is organized as follows. In Sect. 2 we discuss 
the method we use. In Sect. 3 we present the results of the 
application to Planck/LFI. Finally, in Sect. 4 we draw the 
main conclusions of this work. 

2. Method 

2.1. The Statistics of the CMB 

We begin by reviewing some basic notions about the 
statistics of the CMB on the sky sphere. The CMB tem- 
perature field as a function of the direction of observa- 
tion on the sky can be expanded in spherical harmonics 
Ye m (0,d>) as: 



AT(6,ct>) = Y,a im Y 6 



I'm \ 



(1) 



The angular power spectrum of the CMB anisotropy is 
defined as: 



Ce — (|a^m| 2 ) 



(2) 



where the operation (■) represents the average over the 
statistical ensemble. We cannot measure this quantity di- 
rectly, since our own sky is only a particular realization 
of the statistical ensemble. The best possible unbiased es- 
timator of the power spectrum can be constructed by re- 
placing the ensemble average with an average over the 
2£ + 1 independent a,f m coefficients available for each £: 



Cf 



21 



+ i ^ 



(3) 



We have used the superscript S to indicate that this is 
just an estimator of the underlying Ci obtained from our 
particular sky realization. If each ai m coefficient is a zero- 
mean Gaussian variable, as commonly assumed, the Cf 



estimates follows a x 2 distribution with 2£ + 1 degrees of 
freedom (DOF) and rms (e.g. Knox 1995): 



AC? = C, 



21+1 



(4) 



If the CMB fluctuations AT are not measured over 
the entire sky, then it is not possible to use the spherical 
harmonics Yg m as a complete basis to perform the mul- 
tipole expansion. As a consequence, the power spectrum 
estimated using formula (||) will be biased, and values cor- 
responding to different £ will be correlated. The effect of 
incomplete sky coverage on the Ci statistics has been stud- 
ied by Wandelt et al. (2001), Oh et al. (1999), Mortlock et 
al. (2002). If W(6, 4>) = Et <m wi m Y tm {9, 0) is the weight- 
ing function describing a given sky coverage, and: 

m 

then it can be shown (see e.g. Hivon et al. 2001) that: 
(Cf)=^M«,Q (6) 



where: 



and 



21'- 
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(7) 



is the Wigner 3 — j symbol. 





An approximate way to model the effect of partial sky 
coverage is to change the number of DOF of the Ci dis- 
tribution from 2£ + 1 to (2£ + l)/ s ky, where / s ky is the 
observed fraction of the sky. The Cg estimates get biased 
by a / s ky factor, and the rms becomes (Scott et al. 1994; 
Hobson & Magueijo 1996): 



AC; 



(2£+l)/sky 



(8) 



Residual correlations in multipole space can be minimized 
by estimating the average power spectrum in bands of 
appropriate width A£. 

Since the applications shown in this paper will regard 
Cf, estimation from nearly full-sky maps (/sky — 1), an ex- 
act treatment of the partial sky coverage effect is totally 
superfluous, as the approximation described above yields 
comparable accuracy. Indeed, we have verified that this 
simplified approach provides an excellent approximation 
to the exact treatment even when a considerable fraction 
of the sky remains unobserved. We have applied a rather 
severe Galactic cut (20° symmetric around the equator, 
corresponding to / s ky ~ 65.8%) to 100 Monte Carlo simu- 
lations of maps containing only CMB signal, and extracted 
the corresponding Cf. Figure [l] shows the results of this 
test. When the Ci estimates are corrected for the / s k y bias 
they are nearly undistinguishable from the true underlying 
theoretical model. From the same Figure it is also appar- 
ent that Eq. (^]) provides a very good approximation to 
the true dispersion of the estimated Ct. 
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Fig. 1. The effect of incomplete sky coverage on the power 
spectrum estimates. The blue dots are Q// S ky extracted 
from 100 Monte Carlo simulated maps containing only 
CMB signal, after removing a symmetric strip of ±20° 
around the Galactic equator (/ s ky = 65.8%). The black 
dotted curves are the input theoretical Ci used in the sim- 
ulations and the upper and lower la bounds obtained from 
Eq. (0). The blue continuous curves are the mean Cg and 
upper and lower la bounds estimated from the MC simu- 
lations. This shows that Eq. (||) is a good approximation 
even when a substantial fraction of the sky is unobserved. 



We then decided, for the present paper, to neglect the 
complications of an exact treatment of the partial sky cov- 
erage effect. We emphasize, however, that this does not 
limit the applicability of the method, since we are able to 
treat exactly the case of arbitrary sky coverage (for exam- 
ple for very small, irregularly shaped patches, f s k y <C 1, as 
those observed with balloon experiments) by adopting the 
correct formalism (as done, e.g., by Hivon et al. (2001)). 

2.2. Instrumental Noise Contribution to CMB Power 
Spectrum Measurements 

The measurements performed by a given experiment can 
be thought of as a superposition of sky signal (S) and 
instrumental noise (N), which are independent statisti- 
cal processes. When a map is produced from the obser- 
vations, residual noise is left in the pixelized data. This 
noise contribution is minimal when a minimum variance 
map-making (such as the IGLS) is used to obtain the map. 
We can write the observation in pixel p as: 



(9) 



A direct estimation of the power spectrum from such a 
noisy map would yield: 



r <SN 



2f 



— V 

+ 1 ^ 



Expanding this expression, and averaging over the statis- 
tical ensemble, we obtain (due to the independence of the 
sky and noise processes): 



(cf) = (Q s ) + (Q N ). 

Thus, an unbiased estimator of Cf is: 
Cf = Ct»-{C?). 



(11) 



(12) 



Direct extraction of the angular power spectrum from 
a map is a relatively quick task, that can be performed 
in Afp^ 2 log M p operations using fast spherical harmonics 
transforms (Muciaccia et al. 1997; Gorski et al. 1999). This 
facilitates a fast unbiased estimation of the power spec- 
trum through Eq. (|l2|), where the noise bias (Cf) can be 
evaluated using Monte Carlo techniques as discussed in 
the next Section. 

It can be easily shown that the rms uncertainty of this 
power spectrum estimate is: 



AQ E = (Q + Cf) 



{21 + l)/ sky ' 



(13) 



Note that this uncertainty depends on the unknown 
ensemble average value Cg. An approximate way of cal- 
culating ACf is simply to replace Cg with the estimate 
Cf in Eq. (|l3|). A more rigorous strategy is to set frc- 
qucntist confidence intervals by computing the values of 
Ci which are consistent with the estimate Cf at a given 
confidence level. We have verified that the two approaches 
yield nearly identical la error bars for Cf . 

3. Results 

To determine the noise properties in multipole space, we 
first apply the IGLS map-making procedure on the TOD, 
ending up with a minimum variance map of the CMB 
and an estimate of the true noise properties in time do- 
main. We then produce a number of MC realizations of 
time streams containing only instrumental noise, where 
the noise has the statistical properties measured from the 
data. These simulated time streams naturally incorporate 
all information about the experiment observational strat- 
egy. We apply again the IGLS map-making algorithm to 
these mock data sets, to construct a set of noise maps. 
These noise maps have the same statistical properties of 
the noise which contaminates the real map. They can 
then be used to characterize the noise statistical proper- 
ties in multipole space, for example estimating the mean 



(Cf )mC to be used in Eqs. (|12|) and ((131) 



l tm + a fm) 



(10) 



Figure ^| shows the result of applying this procedure to 
simulated observations of Planck/LFI 100 GHz channel. 
The simulations were produced using the standard Planck 
scanning strategy (i.e. a spinning frequency of 1 r.p.m. and 
85° offset angle between the pointing direction and spin 
axis of the telescope, with the latter always on the ecliptic) 
as well as a realistic model of the receiver noise (i.e. 1// 
noise with /knee = 0.1 Hz and Planck sensitivity goal, see 
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Fig. 2. Monte Carlo estimation of the noise properties in 
multipole space for Planck/LFI. The upper plot shows the 
average instrumental noise contribution to the power spec- 
trum, {C^)mc- The blue curve was obtained from 22 re- 
alizations of maps containing only instrumental noise with 
realistic properties (including time correlations); the green 
curve is from 100 realizations of non uniform uncorrelated 
noise; the orange curve is from a subset of 22 realizations 
of the same non uniform uncorrelated noise maps. The 
simulations were performed using the nominal sensitivity 
of all combined 34 receivers at 100 GHz, for 7 months of 
observation. The lower plot shows the noise contribution 
to the error bars, from the same sets of simulations. 



e.g. Natoli et al. (2002)). A comparison was made with the 
case of non uniform uncorrelated noise. It is apparent that 
such a simple model is not accurate enough to characterize 
the noise properties in multipole space, particularly for 
small values of i (£ < 100). It is interesting to observe 
that, since we only need to estimate an average from the 
Monte Carlo, i.e. the Cf to be used in Eqs. ( p"2| ) and 
(|l3|), convergence is achieved after a fairly small number 
of simulations. This is evident from Fig. 0. 

Figure [| shows histograms of the noise angular power 
spectrum Cf , for four multipole values, obtained from 
1 000 Monte Carlo realizations Planck/LFI maps contain- 
ing only instrumental noise. In order to speed up the pro- 
duction of simulated maps, we used an inhomogeneous 
white noise model in place of the correlated noise used 
in the rest of the paper. This should not alter the re- 
sults of this test, since the deviation from the white noise 
behaviour is relevant only at low multipoles (see Fig. ^), 
where the noise contribution is sub-dominant with respect 
to the cosmic variance. The comparison with the theo- 
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Fig. 3. Distribution of the noise contribution in multipole 
space. The plots show the histogram of Cf , normalized to 
its average value, extracted from 1 000 Monte Carlo simu- 
lations Planck/LFI maps, containing only inhomogeneous 
uncorrelated instrumental noise, for four values of I (from 
left to right and from top to bottom: t =50, 200, 400, 
1000). Also shown is the \ 2 distribution with (21 + l)/ s ky 
degrees of freedom (continuous line) and the correspond- 
ing la symmetric error bar, given by y // 2/((2£ + l)/ s ky) 
(blue lines). 



retical distribution, i.e. a x 2 with (2£ + l)/ s ky degrees of 
freedom, shows that use of Eq. ( |l3| ) yields nearly optimal 
error bars for the power spectrum estimates. 

Figures |] and ^| show the results of applying our 
method to extract the CMB power spectrum from a sim- 
ulation of Planck/LFI observations at 30 and 100 GHz. 
In simulating the observations we fully took into account 
the real Planck scanning strategy. The 85° offset between 
the pointing direction and spin axis of the telescope leaves 
unobserved two areas of about 80 square degrees around 
the ecliptic poles. We modeled the optical response of the 
instrument as a symmetric Gaussian beam with FWHM 
of 33 arcminutes for the 30 GHz detectors and of 10 
arcminutes for the 100 GHz detectors. The simulated 
TOD for the 30 GHz channel has Md ~ 10 9 time sam- 
ples (14 months of observation), which are mapped into 
Mp ~ 7 x 10 5 pixels. These numbers turn into Md ~ 2 x 10 9 
(7 months of observation^) and Af p ~ 3 x 10 6 for the 100 
GHz channel (we remind that IGLS map-making scales 



3 Our I/O routine is currently able to manage only 32 bit 
integer. This limits the maximum number of time samples we 



could manage to 2 
being removed. 



2 x 10 9 . This technical limitation is now 



A. Balbi et al.: CMB Power Spectrum Estimation for the Planck Surveyor 



5 



approximately linearly with A/d). Going from the simu- 
lated data to the power spectrum estimates took about 
12 h for the 30 GHz channel, and about 18 h for the 100 
GHz channel, using a parallel implementation of the IGLS 
algorithm running on 16 processors of the Origin 3000 su- 
percomputer at Cineca. These times are dominated by the 
production of the MC noise maps, 22 in our case. This 
number is enough to produce good estimates for the av- 
erage noise contribution (see again Fig. |^) . For a detailed 
discussion about memory requirements, CPU timing, and 
code scalability of our implementation of the IGLS map- 
making algorithm, see the paper by Natoli et al. (2001). 

Comparison of our Ci estimates with the theoretical 
input model yields a x 2 /d.o.f. of 1.06 and 1.10 for the 30 
and 100 GHz case, respectively. As mentioned in Section 



2.1 



we also estimated the power C& = + l)Ce/2ir 

in bands b of width A£ — 20, in order to reduce small 
spurious correlation in multipole space. When we use a 
simple cubic spline algorithm to interpolate these band- 
power estimates, the resulting curve is nearly undistin- 
guishable from the theoretical input model (see bottom 
panel in Figs. [| and |^). 



4. Conclusions 

We have shown that IGLS map-making can be successfully 
applied, in conjunction with MC techniques, to the prob- 
lem of estimating the angular power spectrum from CMB 
maps. The method discussed in this work is fast enough 
to be already applicable to megapixel maps such as those 
expected from the Planck Surveyor. No manipulation of 
the time stream (i.e. high-pass filtering) is required by 
this method. Furthermore, no unrealistic simplification of 
the instrumental noise behavior is needed, since the noise 
properties are estimated directly from the data. Our es- 
timated noise angular power spectrum contains the right 
information on noise correlation, as well as on the details 
of scanning strategy, even if we never use explicitly the 
full pixel covariance matrix. Although we only presented 
in this paper an application to the case of nearly full-sky 
maps from the Planck Surveyor, the approach described 
here can be used for any arbitrary scanning strategy and 
sky coverage (for example to analyze balloon data), as 
long as the spectral information on the spatial window of 
the observation is taken into account. 

Finally, we would like to mention a few possible de- 
velopments of the work described in this paper. We con- 
sidered a symmetric beam approximation in our analysis. 
This approximation may prove to be too simplistic when 
dealing with real experiments. Thus, we are currently gen- 
eralizing our map-making algorithm to deal with asym- 
metric beam patterns. We also point out that this power 
spectrum estimation technique is easily generalizable to 
CMB polarization maps, since an implementation of the 
IGLS map-making algorithm for polarization observations 
exists (see Balbi et al. 2002). We shall discuss this topic 
in a forthcoming paper. 
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Fig. 4. CMB Power spectrum estimation for Planck/LFI 
30 GHz channel. The simulated data were produced as- 
suming the nominal sensitivity from all 4 combined ra- 
diometers and 14 months of observation (full mission). 
The top panel shows the power spectrum estimated di- 
rectly from the map {Cf N , green points), the MC esti- 
mation of the noise contribution ((C^)mCi blue curve), 
and the recovered signal power spectrum (Cf , red points). 
The bottom panel shows the estimated band-powers Cf, = 
J2ieb^(^ + l)Ce/2n and their la error bars (red boxes). 
The orange curve is the result of a cubic spline interpola- 
tion of the band-power estimates. In both panels the black 
curve is the input theoretical model used in the simulation 
(a COBE- normalized standard CDM). 
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Fig. 5. Same as in Fig. |J, but for the 34 combined ra- 
diometers of Planck/LFI 100 GHz channel and 7 months 
of observation (half mission) . 
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